m6A mRNA Methylation Analysis Provides Novel Insights into Pigmentation in Sheep Skin

ABSTRACT N6-methyladenosine (m6A) is the most universal post-transcriptional modification of mRNA which may play important roles in verious species. However, the potential roles of m6A in the pigmentation of skin are not completely understood. To explore the role of m6A modification in pigmentation of sheep skin, we used MeRIP-seq and RNA-seq to profile the skin transcriptome in black and white coat color (n=3). Our results showed that an average of 7701 m6A peaks were obtained for all samples and the average length was 305.89 bp. The GGACUU sequence was the most enrichment motif and shared in black skin and white skin. The m6A peaks were mainly enriched in the CDS, 3'UTR and 5'UTR, especially in CDS region near the stop codon of the transcript. 235 signiﬁcantly differential peaks were found in black skin vs. white skin. The KEGG signaling pathways of downregulated and upregulated m6A peaks were mainly enriched in AGE-RAGE signaling pathway in diabetic complications, Viral carcinogenesis, Transcriptional misregulation in cancer, ABC transporters, Basal transcription factors and Thyroid hormone synthesis (P value <0.05). For RNA-seq, 71 differently expressed genes (DEGs) were scanned in black skin vs. white skin. DEGs were significantly enriched in tyrosine metabolism, melanogenesis, neuroactive ligand-receptor interaction pathway (P value <0.05). Combined m6A-seq and RNA-seq analysis showed that the hyper-up genes and hypo-up genes were both enriched in ErbB signaling pathway (P value <0.05). In conclusion, it provide a basis for further research into the functions of m6A methylation modiﬁcations in pigmentation.


Introduction
Coat colour is a direct reflection of pigmentation in fleece-producing animals. In mammals, the main pigments of the skin and hair are melanins, which are synthesized in melanosomes of melanocytes and then secreted into keratinocytes [1]. Multiple genes and signalling pathways involves in melanin synthesis in mammals, such as tyrosinase (encoded by TYR) catalyse the oxidation of tyrosine (or of L-3,4-dihydroxyphenyla-lanine [L-dopa]) to dopaquinone, which is the initial reaction in melanin synthesis [2]. TYRP1 (TYRrelated protein 1) and DCT/TYRP2 (dopachrome tautomerase) modulate the eumelanin synthesis [3,4]. Pre-melanosome protein (encoded by PMEL) is the component of the intraluminal fibrous sheets on which eumelanins are deposited in melanosomes [5]. SLC45A2 directly impact the pH of maturing melanosomes [6].
To explore the role of m 6 A modification in pigmentation of sheep skin, we performed N6methyladenosine sequencing (m 6 A-seq) and RNA sequencing (RNA-seq) to investigate differentially methylated genes (DMGs) and differentially expressed genes (DEGs) in black and white skin of sheep. Our results provide a theoretical basis for further research into the molecular mechanisms of m 6 A modification in pigmentation.

Materials and methods
All related experiments involving sheep were conducted in strict compliance with relevant guidelines set by the Ethics Committee of Tongren University, China (Approval ID: TREDU2022-016).

Animals and sample collection
Small Tailed Han Sheep used in this study were raised under the same conditions at Taigu Haihong Animal Husbandry Co. LTD. (Taigu, China). Three one-year-old sheep similar in size with black-and-white coat colour were selected. The black hair and adjacent white hair were cut off, and the black skin and adjacent white skin were taken using skin biopsy borer with a diameter of 1 cm. Four skin pieces of each colour were taken from each sheep. All samples were put into 1.5 mL centrifuge tubes, labelled and stored in liquid nitrogen for RNA extraction.

RNA isolation and fragmentation
Total RNA from 6 samples (1 black and 1 white skin tissue each sheep) was extracted using TRlzol™ Reagent (Thermo Fisher Scientifi, MA, USA) according to the manufacturer's instructions and genomic DNA was removed with DNase I (Roche Diagnostics, IN, USA). A260/A280 > 1.9, A260/A230 > 1.7 was used to evaluate RNA purity using a NanoDrop ND-1000 instrument (NanoDrop, DE, USA). The concentration of total RNA was measured by Qubit RNA HS Assay Kit (Thermo Fisher Scientific, MA, USA). The purified mRNA (25 μg) was randomly broken into ~150 nt fragments using RNA Fragmentation Buffer (100 mM TrisHCl, 100 mM ZnCl 2 ) at 70°C for 6 min. 98% of the fragmented mRNA was used for immunoprecipitation (IP), and the rest (2%) for IP control (Input).

RNA immunoprecipitation
m 6 A MeRIP is based on the previously described m 6 A-seq protocol [21] with several modifications: 30 μL of protein A magnetic beads (Thermo Fisher Scientific, MA, USA) were washed twice by IP buffer (150 mM NaCl, 10 mM Tris-HCl [pH 7.5], 0.1% IGEPAL CA-630), resuspended in 500 μL of IP buffer, and tumbled with 5 μg anti-m 6 A antibody (Synaptic Systems, Gottingen, Germany) at 4°C for 6 h. It was then incubated in 500 μL IP buffer with 5 μL of RNasin Plus RNase Inhibitor (40 U/μL, Promega, WI, USA) at 4°C for 2 h. After extensive washing, the m 6 A-enriched fragmented RNA was eluted by 200 μL of RLT buffer supplied in RNeasy Mini Kit (QIAGEN, Dusseldorf, Germany) for 2 min at room temperature. Thereafter, supernatant was collected and 400 μL of 100% ethanol was added, then m 6 A-enriched RNA was purified using an RNeasy MiniElute spin column (QIAGEN, Dusseldorf, Germany). Finally, The m 6 A-enriched RNA was eluted with 14 μL ultrapure H 2 O.

M 6 A library preparation and sequencing
The amount of 2 μL eluted RNA and input RNA was reverse transcribed with High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific, MA, USA). Transcriptome-wide interrogation was pursued by deep sequencing using SMARTer Stranded Total RNA-Seq Kit version 2 (Pico Input Mammalian, Takara/Clontech, Osaka, Japan) on Illumina Novaseq 6000 platform. 150 bp paired-end reads were generated.

RNA-seq
A total amount of 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample.

Bioinformatic analysis of m 6 A-seq
Trimmomatic software (v0.32) were used to remove adapter and low quality reads [22]. Quality distribution plot and base content distribution were generated by FastQC (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/). The software STAR was used to align reads to the sheep reference genome (Oar_v4.0, https://www.ncbi.nlm.nih.gov/assembly/GCF_ 000298735.2/). Uniquely mapped reads were used for subsequent analysis. The R software package MetPeak (Cutoff threshold: PEAK_CUTOFF_P value = 0.05, FOLD_ENRICHMENT = 1) was used to call peak and IGV software (http://www.igv.org) were used to visualize. The R software package Guitar was used to calculate the densty plot of peaks in 3'UTR, CDS and 5'UTR. ChIPseeker (https://bioconductor.org/ packages/ChIPseeker) and HOMER (http://homer. ucsd.edu/homer/motif) softwares were used to annotate the peaks and perform motif analysis, respectively. The MeTDiff software was used to analyse the difference the MeRIP-seq data between black and white skin (P value < 0.05 and Log 2 FC >1 as upregulated peak, P value < 0.05 and Log2 FC <-1 as downregulated peak). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of m 6 A modified genes (The 3'UTR of differentially methylated mRNA was modified by m 6 A) were performed by the DAVID database (http://david.abcc. ncifcrf.gov/) and the KEGG database (https://www. kegg.jp/), respectively.

Bioinformatic analysis of RNA-seq
Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts to obtain clean data (clean reads), which were used for further analysis. STAR was used to align clean reads to reference genome. HTSeq v0.6.0 was used to count the reads numbers mapped to each gene. And then FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs) was calculated based on the length and reads count of genes. The differentially expressed genes (DEGs) were filtered by DESeq2 algorithm with the criteria of |log 2 FC| > 1 and P value < 0.05. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of differentially expressed genes were performed by the DAVID database (http:// david.abcc.ncifcrf.gov/) and the KEGG database (https://www.kegg.jp/), respectively. Fisher's exact test was applied to identify the significant GO categories, significant pathway and FDR was used to correct the P values.

Real-Time quantitative PCR (RT-Qpcr)) and m 6 A IP (MeRIP) followed RT-Qpcr (MeRIP-Qpcr)
IP RNA and Input RNA (1 μg) from each sample were synthesized cDNA by PrimeScript™ RT Master Mix (Perfect Real Time) (TAKARA, Dalian, China). Thereafter, a TB Green® Fast qPCR Mix (TAKARA, Dalian, China) and a LightCycler 480II (Roche, Basel, Switzerland) were used to perform RT-qPCR and MeRIP-qPCR. The primers are presented in Table 1, β-actin served as internal control. The relative expression of genes were determined by the 2 -ΔΔCt .

MeRIP-seq summary
To explore the role of m 6 A modification in pigmentation of sheep skin, the black and white skins of three sheep were collected for m 6 A sequencing. The IP libraries (MeRIP-seq) and the input libraries (RNA-seq) of black skin and white skin were constructed. As a result, a total of 56,239,292-58,047,254 and 61,572,416-65,348,920 raw reads in the IP libraries and in the input libraries were obtained for black skin, respectively. Meanwhile, a total of 53,302,052-59,653,724 and 55,037,238-64,703,960 raw reads in the IP libraries and in the input libraries were obtained for white skin, respectively. An average of more than 50,143,426 clean reads per sample were obtained. The Q20 and Q30 values were at least 97.15% and 91.74% for all sample, respectively ( Table 2). In addition, more than 50% clean reads were aligned to the sheep reference genome. At least 18,019,944 clean reads were uniquely mapped and the percent of uniquely mapped reads >70% for all samples (Supplementary Table S1).

General features of sheep m 6 A methylation
An average of 7701 peaks were obtained for all samples and the average length was 305.89 bp (Supplementary Table S2, S3). After merging the three replicates, 9354 m 6 A modified peaks and 5217 genes were found in black skin. Meanwhile, 9213 m 6 A modified peaks and 5253 genes were detected in white skin (Figures 1a, 1b). The count of m 6 A peaks in one modified gene was in the range of 1 to 30, and the average was 1.2 to 1.32 m 6 A peaks in six samples. Meanwhile, 80.70% of the modified genes had only one or two m 6 A peaks, the remaining genes (19.30%) contained three or more peaks in sheep skin (Figure 1c). Furthermore, we analysed the distributions of peaks on the sheep chromosomes in skin. Interestingly, there were more peaks on chromosomes 1, 2, 3, 5, 11, and 14 than on other chromosomes, and the most peaks were distributed on chromosome 1 in black skin and chromosome 3 in white skin (Figure 1d,e). To identify RRACH (R, purine; A, m 6 A; H, nonguanine base) motif of sheep skin, HOMER software was used. The results showed that more than 15 motifs are identified in each sample and the five motifs with the smallest P value from each group were used for subsequent analysis (Figure 1f). The GGACUU sequence was shared in black and white skin, which inversely complemented two miRNA (hsa-miR-302e and hsa-miR-2114) seed sequences.

Topological pattern of sheep m 6 A methylation
The distributions of peaks in gene functional elements were analysed. As a result, the m 6 A peaks in black skin and white skin were mainly enriched in the CDS, 3'UTR and 5'UTR, especially in CDS region near the stop codon of the transcript ( Figure 2a). To further evaluate the distribution of m 6 A modified peak, the transcript was divided into the 5'UTR, 3'UTR, 1st exon and other exon. As a result, the highest enrichment of m 6 A modified peaks was in other exon, followed by 3'UTR, 1st exon and 5'UTR ( Figure 2b).

Differential m 6 A peaks between black and white skin
MeTDiff software was used to analyse differential m 6 A peaks between black skin and white skin based on P value < 0.05 and | Log2 FC | >1. The result showed that 235 differential peaks and 226 genes (differentially methylated genes, DMGs) were scanned, of which 134 significantly upregulated peaks in 127 genes and 101 significantly downregulated peaks in 99 genes were found in black skin vs. white skin (Figure 3a and Supplementary Table S4). The top 10 upregulated genes and downregulated genes are present ( Table 3).
The distribution of differential m 6 A peaks on chromosomes showed that the downregulated peaks were most enriched on chromosome 2 (11 differential peaks), while the upregulated peaks were most enriched on chromosome 1 (18 differential peaks). The number of the downregulated peaks and upregulated peaks on chromosome 1 were 6 and 18 was most different ( Figure 3B). The distribution of differential peaks in genes was counted. Only two genes (CSPG4, MAML3) contained two downregulated peaks, but five genes (ADAMTS1, COL6A2, GAB2, MAP1B, TNRC18) contained two upregulated peaks. Even PLEC gene owned 3 upregulated peaks. The remaining genes (97.98% downregulated genes and 95.28% upregulated genes) all have a single differential peak (Supplementary Table S4).

GO and KEGG analysis of genes presenting differential m 6 A Peaks
To explore the role of m 6 A modification in sheep skin and its relationship with pigmentation, the functions of m 6 A modified genes were analysed using the DAVID database and the KEGG database. According to GO terms, the functions were divided into three categories: the biological process (BP), cellular component (CC), and molecular function (MF) categories. Downregulated m 6 A peaks were significantly enriched in 97 BP, 18 CC and 21 MF GO terms, and Upregulated m 6 A peaks were significantly enriched in 21 BP, 7 CC and 8 MF GO terms, the top 15 of three categories were showed in Figure 4. GO terms related to pigmentation involve in regulation of synaptic transmission, dopaminergic, melanosome assembly, regulation of dopamine metabolic process, melanosome organization, negative regulation of kinase activity and melanosome membrane. The KEGG signalling pathways of downregulated m 6 A peaks were mainly enriched in AGE-RAGE signalling pathway in diabetic complications, Viral carcinogenesis, Transcriptional misregulation in cancer, ABC transporters, Basal transcription factors (P value < 0.05). Meanwhile,   6 A peaks were mainly enriched in Thyroid hormone synthesis (P value < 0.05) ( Figure 5).

Differently expressed genes (RNA-seq) in black and white skin
The RNA-seq data for each sample were used to analyse differently expressed genes between black skin and white skin. A total of 71 DEGs were scanned, among which 27 DEGs were downregulated and 44 DEGs were upregulated in black skin vs. white skin (Figure 6a, Supplementary Table  S5). The top 10 upregulated genes and downregulated genes are listed in Table 4. The functions of these DEGs were analysed and the GO enrichment and KEGG pathway were evaluated. As a result, the top 15 GO terms of biological process (BP), cellular component (CC) and molecular function (MF) were exhibited in Figure 6b. The DEGs were significantly enriched in melanin biosynthetic process, response to blue light, transmembrane transport, melanosome organization, positive regulation of protein kinase A signalling, developmental pigmentation, melanosome and L-dopa decarboxylase activity, which were closely related to pigmentation of skin. KEGG analysis showed that DEGs were significantly enriched in Tyrosine metabolism (4 DEGs), Melanogenesis (4 DEGs), Neuroactive ligand-receptor interaction (6 DEGs

Combined m 6 A-seq and RNA-seq analysis
To further determine the effect of m 6 A modification on pigmentation in sheep skin, the relationship of m 6 A modified genes and DEGs were analysed according to m 6 A-seq data and RNA-seq data. As a result, BIIIB4 was the unique overlap gene which was downregulated in both m 6 A methylation and mRNA expression level in black skin vs white skin in one of the three sheep (Figure 7a). In addition, With |log 2 FC|>0 as the dividing line of m 6 A methylation and |log 2 FC|>1 as the threshold of mRNA differential expression, a total of 27 genes were divided into four groups: 9 hypermethylated and downregulated genes (hyper-down genes), 6 hypermethylated and upregulated genes (hyper-up genes), 11 hypomethylated and downregulated genes (hypodown genes), 1 hypomethylated and upregulated genes (hypo-up genes) (Figure 7b and Supplementary Table S6). All genes were used for function analysis using KEGG database, which revealed that the hyper-up genes were mainly enriched in Gap junction, ErbB signalling pathway, Inflammatory mediator regulation of TRP channels and Serotonergic synapse (P value < 0.05) (Figure 7c). In addition, the hypo-up genes were mainly related to ErbB signalling pathway and Amyotrophic lateral sclerosis (ALS) (P value < 0.05) (Figure 7d). However, the hyper-down genes and hypo-down genes were not significantly enriched in any pathway (P value > 0.05)

Validation of DEGs and DMGs by qPCR and MeRIP-Qpcr
To verify the MeRIP-Seq and RNA-seq sequencing data, we selected six DMGs (CSPG4, MAML2, DENND2B, ADAMTS1, MAP1B, GAB2) and six DEGs (TYRP1, SLC45A2, TYR, TCL1A, ZAR1, MCP-3) randomly and detected their expressions by MeRIP-qPCR and RT-qPCR. The result showed that the expression trends of these genes were consistent with the RNA-seq and m 6 A-seq (Figure 8a,b), which confirmed the accuracy of m 6 A-seq and RNA-seq experiment.

Discussion
Coat colour is an important trait in sheep. Coat colour is determined by the content and composition of melanin. Melanin produced in melanocytes can be divided into two types, eumelanin, which appears black to brown and pheomelanin which appears yellow to reddish brown. The genetic basis and many genes controlling coat colour of sheep has been found, such as melanocortin-1 receptor (MC1R), agouti signalling protein (ASIP), TYR, TYRP1 and microphthalmia transcription factor (MITF) [23].
M 6 A modifications, regulating the stability, splicing, translation, and degradation of mRNAs, may play important roles in growth, reproduction, nerve development, fat metabolism, immune responses, tumour invasion and other physiological processes [15,16]. However, the function of m 6 A modifications on pigmentation is unclear. In this study, the m 6 A modifications of black skin and white skin of sheep were detected, 7358 peaks that were in common in the black skin and white skin, which accounted for 79.9% of all the peaks. This means that most m 6 A modifications are designed to maintain cellular metabolism. The average number of m 6 A peaks per gene was 1.2 to 1.32 in six samples, which is similar to that in pig liver (1.33-1.42) and mouse liver (1.34) [24,25]. The consensus sequences 'RRACH' were conserved in various species [26], which was also verified in sheep. The GGACUU sequence was the most enriched m 6 A motif in all samples. The distributions of peaks in gene functional elements were analysed. As a result, the m 6 A peaks in black skin and white skin of sheep were mainly enriched in the CDS, 3'UTR and 5'UTR, especially in CDS region near the stop codon of the transcript. The result agrees with that in human, mouse, pig and sheep liver tissue [14,15,25]. The stability, localization, expression, and translation of mRNA were regulated by 3'UTR where multiple RNA-binding proteins bind to plays a regulatory role and protein interaction [27]. Differential m 6 A peaks between black skin and white skin were analysed. The result showed that 235 differential peaks and 226 genes were scanned, of which 134 significantly upregulated peaks in 127 genes and 101 significantly downregulated peaks in 99 genes were found in black skin vs. white skin. Eight DMGs contains two or more peaks. Chondroitin sulphate proteoglycan 4 (CSPG4) is essential to the survival and growth of melanoma tumours by enhancing growth factor receptor-regulated pathways, such as sustained activation of ERK 1,2, and modulating integrin function [28]. ERK pathway was the key intracellular signalling pathway, which was involved in melanogenesis by the activation of signal transduction [29]. ADAMTS1(A Disintegrin-like And Metalloprotease domain with ThromboSpondin type I motifs) is upregulated in black skin compared with white skin is in accord with previous study. ADAMTS genes (ADAMTS1, ADAMTS6 and ADAMTS9) may associate with age-related macular degeneration and MAP1B was higher expression in the macula of human eye [30,31]. GAB2 amplification is critical for melanomas arising from sunprotected sites and in mast cell development GAB2 is required for KitL/c-Kit signalling which is an important signal in melanogenesis [32][33][34]. Mastermind-like 3 (MAML3) are potential therapeutic targets for small cell lung cancer and pancreatic cancer [35]. Collagen VI (COL6A1, COL6A2 and COL6A3) mutations result in disorders abnormal skin findings [36]. But there is no evidence that MAML3, COL6A2, TNRC18, PLEC are related to pigmentation. This study provides a new direction for these genes. Transcriptome of each sample were detected and a total of 71 DEGs were scanned, among which 27 DEGs were downregulated and 44 DEGs were upregulated in black skin vs. white skin. Upregulated genes TYR, TSPAN10, TRPM1, MLANA, KCNJ13, TYRP1, DCT, SLC45A2, SLC24A5, MC1R and PMEL have been proven to participate in pigmentation [37][38][39][40][41][42].
DEGs (RNA-seq) in black skin vs. white skin were enriched GO terms of melanin biosynthetic process, response to blue light, transmembrane transport, melanosome organization, positive regulation of protein kinase A signalling, developmental pigmentation, melanosome and L-dopa decarboxylase activity. KEGG analysis showed that DEGs were significantly enriched in tyrosine metabolism (4 DEGs), melanogenesis (4 DEGs), neuroactive ligand-receptor interaction (6 DEGs). The result is similar to the previous study [41,42]. Meanwhile, GO terms of m 6 A modified genes also involved in pigmentation, such as regulation of synaptic transmission, dopaminergic, melanosome assembly, regulation of dopamine metabolic process, melanosome organization, negative regulation of kinase activity and melanosome membrane, which is similar to GO enrichment of DEGs. Combined analysis of m 6 A-seq and RNA-seq showed that both the hyper-up genes and hypo-up genes enriched in ErbB signalling pathway which are required to promote normal pigment cell and pigment pattern development in vivo [43].

Conclusion
In this study, we detected how m 6 A methylation is modified in black skin and white skin of sheep, and m 6 A modification may play an important role in pigmentation of skin in sheep by participating in AGE-RAGE signalling pathway, ABC transporters, Basal transcription factors and Thyroid hormone synthesis. It provides a basis for further research into the functions of m 6 A methylation modifications in pigmentation.